Load the necessary libraries
library(car) # for regression diagnostics
library(broom) # for tidy output
library(ggfortify) # for model diagnostics
library(sjPlot) # for outputs
library(knitr) # for kable
library(effects) # for partial effects plots
library(ggeffects) # for partial effects plots
library(emmeans) # for estimating marginal means
library(MASS) # for glm.nb
library(MuMIn) # for AICc
library(tidyverse) # for data wrangling
library(nlme)
library(lme4)
library(glmmTMB)
library(broom.mixed)
library(glmmTMB) # for glmmTMB
library(DHARMa) # for residuals and diagnostics
library(performance) # for diagnostic plots
library(see) # for diagnostic plots
library(patchwork) # for multiplot layouts
To investigate differential metabolic plasticity in barramundi (Lates calcarifer), Norin, Malte, and Clark (2015) exposed juvenile barramundi to various environmental changes (increased temperature, decreased salinity and increased hypoxia) as well as control conditions. Metabolic plasticity was calculated as the percentage difference in standard metabolic rate between the various treatment conditions and the standard metabolic rate under control conditions. They were interested in whether there was a relationship between metabolic plasticity and typical (control) metabolism and how the different treatment conditions impact on this relationship.
A total of 60 barramundi juveniles were subject to each of the three conditions (high temperature, low salinity and hypoxia) in addition to control conditions. Fish mass was also recorded as a covariate as this is known to influence metabolic parameters.
Barramundi
Sampling design
Format of norin.csv data files
| FISHID | MASS | TRIAL | SMR_contr | CHANGE |
|---|---|---|---|---|
| 1 | 35.69 | LowSalinity | 5.85 | -31.92 |
| 2 | 33.84 | LowSalinity | 6.53 | 2.52 |
| 3 | 37.78 | LowSalinity | 5.66 | -6.28 |
| .. | .. | .. | .. | .. |
| 1 | 36.80 | HighTemperature | 5.85 | 18.32 |
| 2 | 34.98 | HighTemperature | 6.53 | 19.06 |
| 3 | 38.38 | HighTemperature | 5.66 | 19.03 |
| .. | .. | .. | .. | .. |
| 1 | 45.06 | Hypoxia | 5.85 | -18.61 |
| 2 | 43.51 | Hypoxia | 6.53 | -5.37 |
| 3 | 45.11 | Hypoxia | 5.66 | -13.95 |
| FISHID | Categorical listing of the individual fish that are repeatedly sampled |
| MASS | Mass (g) of barramundi. Covariate in analysis |
| TRIAL | Categorical listing of the trial (LowSalinity: 10ppt salinity; HighTemperature: 35 degrees; Hypoxia: 45% air-sat. oxygen. |
| SMR_contr | Standard metabolic rate (mg/h/39.4 g of fish) under control trial conditions (35 ppt salinity, 29 degrees, normoxia) |
| CHANGE | Percentage difference in Standard metabolic rate (mg/h/39.4 g of fish) between Trial conditions and control adjusted for 'regression to the mean'. |
norin <- read_csv("../public/data/norin.csv", trim_ws = TRUE)
glimpse(norin)
## Rows: 180
## Columns: 5
## $ FISHID <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
## $ MASS <dbl> 35.69, 33.84, 37.78, 26.58, 37.62, 37.68, 30.62, 50.37, 24.9…
## $ TRIAL <chr> "LowSalinity", "LowSalinity", "LowSalinity", "LowSalinity", …
## $ SMR_contr <dbl> 5.847466, 6.530707, 5.659556, 6.278200, 4.407336, 4.818589, …
## $ CHANGE <dbl> -31.919389, 2.520929, -6.284968, -4.346675, -3.071329, -15.0…
Since the response variable is a difference between the SMR under control conditions compared to the specific trial conditions (either Low salinity, High temperature or Hypoxia), the only possible distribution we can use is a Gaussian. Well to clarify this, we could use a t-distribution, but this is not really available for frequentist mixed effects routines in R.
Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]
where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of temperature and (centred) mean fish size on SDA peak. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual fish.
Lets begin by declaring the categorical predictors and random effects as factors.
norin <- norin %>% mutate(
FISHID = factor(FISHID),
TRIAL = factor(TRIAL)
)
Each of the fish were exposed to each of the three trials (treatments). So lets explore the distributions of responses within each of these trials.
ggplot(norin, aes(y = CHANGE, x = TRIAL)) +
geom_boxplot()
Conclusions:
The researchers considered that physiological plasticity might be effected by the fishes basal metabolic rate. For example, a fish with relatively high metabolism might have less scope to increase this metabolism than a fish with a lower metabolism. Therefore, the SMR under control conditions (just prior to the specific trial) could be added as a continuous covariate.
Doing so would introduce two additional assumptions:
ggplot(norin, aes(y = CHANGE, x = SMR_contr, shape = TRIAL, color = TRIAL)) +
geom_smooth(method = "lm") +
geom_point()
ggplot(norin, aes(y = CHANGE, x = SMR_contr, shape = TRIAL, color = TRIAL)) +
geom_smooth() +
geom_point()
It might also be worth exploring how consistent the trial effect is across the different fish. This can give us an idea of whether the addition of a random intercept/slope model would be useful.
ggplot(norin, aes(y = CHANGE, x = as.numeric(FISHID), color = TRIAL)) +
geom_point() +
geom_line()
Conclusions:
Finally, as an acknowledgement that the metabolic response might be influenced
by the mass of the fish, the researchers contemplated including fish Mass as a continuous covariate. There are numerous alternative ways to incorporate covariates that are known to impact on a response.
Lets explore the relationship between the response and Mass separately for each Trial.
ggplot(norin, aes(y = CHANGE, x = MASS, color = TRIAL)) +
geom_point() +
geom_smooth(method = "lm")
Conclusions:
For each of the following modelling alternatives, we will:
## Start with random intercept model
norin.lme1 <- lme(CHANGE ~ 1 + offset(MASS),
random = ~ 1 | FISHID,
data = norin,
method = "REML"
)
## Then random intercept and slope model
norin.lme2 <- lme(CHANGE ~ 1 + offset(MASS),
random = ~ TRIAL | FISHID,
data = norin,
method = "REML"
)
AICc(norin.lme1, norin.lme2)
anova(norin.lme1, norin.lme2)
Conclusions:
## Compare models that estimate partial slope for MASS vs an offset for MASS
## must use ML to compare models that vary in fixed effects
norin.lme2a <- update(norin.lme2,
. ~ TRIAL * scale(SMR_contr, scale = FALSE) + offset(MASS),
method = "ML"
)
## The following does not converge
## norin.lme2b <- update(norin.lme2,
## .~TRIAL*scale(SMR_contr, scale=FALSE) + MASS,
## method = 'ML')
## Use BFGS from optim instead
norin.lme2b <- update(norin.lme2,
. ~ TRIAL * scale(SMR_contr, scale = FALSE) + MASS,
method = "ML",
control = lmeControl(
opt = "optim",
optimMethod = "BFGS"
)
)
norin.lme2c <- update(norin.lme2,
. ~ TRIAL * scale(SMR_contr, scale = FALSE),
method = "ML"
)
AICc(norin.lme2a, norin.lme2b, norin.lme2c)
anova(norin.lme2a, norin.lme2b, norin.lme2c)
## Now that we have decided on the structure, we need to ensure that the
## model has been fit with REML
norin.lme2c <- update(norin.lme2c, method = "REML")
Conclusions:
## Start with random intercept model
norin.lmer1 <- lmer(CHANGE ~ 1 + offset(MASS) + (1 | FISHID),
data = norin,
REML = TRUE
)
## Then random intercept and slope model
norin.lmer2 <- lmer(CHANGE ~ 1 + offset(MASS) + (TRIAL | FISHID),
data = norin,
REML = TRUE,
control = lmerControl(check.nobs.vs.nRE = "ignore")
)
AICc(norin.lmer1, norin.lmer2)
anova(norin.lmer1, norin.lmer2)
Conclusions:
## Compare models that estimate partial slope for MASS vs an offset for MASS
## must use ML to compare models that vary in fixed effects
norin.lmer2a <- update(norin.lmer2,
. ~ TRIAL * scale(SMR_contr, scale = FALSE) + offset(MASS)
+ (TRIAL | FISHID),
data = norin,
REML = FALSE
)
norin.lmer2b <- update(norin.lmer2,
. ~ TRIAL * scale(SMR_contr, scale = FALSE) + MASS
+ (TRIAL | FISHID),
data = norin,
REML = FALSE
)
norin.lmer2c <- update(norin.lmer2,
. ~ TRIAL * scale(SMR_contr, scale = FALSE)
+ (TRIAL | FISHID),
data = norin,
REML = FALSE
)
## Compare these models via AICc
AICc(norin.lmer2a, norin.lmer2b, norin.lmer2c)
## Alternatively, we can use sequential Likelihood Ratio Tests (LRT)
anova(norin.lmer2a, norin.lmer2b, norin.lmer2c)
## Now that we have decided on the structure, we need to ensure that the
## model has been fit with REML
norin.lmer2c <- update(norin.lmer2c, REML = TRUE)
Conclusions:
## Start with random intercept model
norin.glmmTMB1 <- glmmTMB(CHANGE ~ 1 + offset(MASS) + (1 | FISHID),
data = norin,
REML = TRUE
)
## Then random intercept and slope model
norin.glmmTMB2 <- glmmTMB(CHANGE ~ 1 + offset(MASS) + (TRIAL | FISHID),
data = norin,
REML = TRUE
)
## norin.glmmTMB2 <- glmmTMB(CHANGE ~ 1 + offset(MASS) + (TRIAL|FISHID),
## data = as.data.frame(norin),
## REML = TRUE,
## control=glmmTMBControl(optimizer=optim,
## optArgs=list(method='BFGS')
## )
## )
norin.glmmTMB2 <- glmmTMB(CHANGE ~ 1 + offset(scale(MASS, scale = FALSE)) + (TRIAL | FISHID),
data = norin,
REML = TRUE
)
AICc(norin.glmmTMB1, norin.glmmTMB2)
anova(norin.glmmTMB1, norin.glmmTMB2)
Conclusions:
## Compare models that estimate partial slope for MASS vs an offset for MASS
norin.glmmTMB2a <- glmmTMB(CHANGE ~ TRIAL * scale(SMR_contr, scale = FALSE) +
offset(MASS) +
(TRIAL | FISHID),
data = norin,
REML = FALSE
)
## OR with update
norin.glmmTMB2a <- update(norin.glmmTMB2,
. ~ TRIAL * scale(SMR_contr, scale = FALSE) +
offset(scale(MASS, scale = FALSE)) +
(TRIAL | FISHID),
REML = FALSE
)
## now with MASS as a partial slope
norin.glmmTMB2b <- glmmTMB(CHANGE ~ TRIAL * scale(SMR_contr, scale = FALSE) +
scale(MASS, scale = FALSE) +
(TRIAL | FISHID),
data = norin,
REML = FALSE
)
## OR with update
norin.glmmTMB2b <- update(norin.glmmTMB2,
. ~ TRIAL * scale(SMR_contr, scale = FALSE) +
scale(MASS, scale = FALSE) +
(TRIAL | FISHID),
REML = FALSE
)
## Finally without MASS
norin.glmmTMB2c <- glmmTMB(CHANGE ~
TRIAL * scale(SMR_contr, scale = FALSE) +
(TRIAL | FISHID),
data = norin,
REML = FALSE
)
## I have explored a range of alternative optimizers etc and
## none of them resolve the issue.
## The current model will be good enough to explore AICc and it turns out that the model
## converges fine once we run with REML
## Note the following does resolve it (not sure why), yet causes downstream issues
## norin.glmmTMB2c <- glmmTMB(CHANGE~
## TRIAL*scale(SMR_contr, scale=FALSE) +
## (TRIAL|FISHID),
## data = norin,
## control=glmmTMBControl(profile=TRUE),
## REML=FALSE)
## ## OR with update
## norin.glmmTMB2c <- update(norin.glmmTMB2,
## .~TRIAL*scale(SMR_contr, scale=FALSE) +
## (TRIAL|FISHID),
## control=glmmTMBControl(profile=TRUE),
## REML=FALSE)
## Compare these models via AICc
AICc(norin.glmmTMB2a, norin.glmmTMB2b, norin.glmmTMB2c)
## Now that we have decided on the structure, we need to ensure that the
## model has been fit with REML
## It turns out that when the model withough MASS is run using REML, the
## response type residuals are tiny - not sure why.
## We will instead go with the offset(MASS) model (2a)
norin.glmmTMB2c <- update(norin.glmmTMB2a, REML = TRUE)
Conclusions:
Conclusions:
plot_model(norin.lme2c, type = "diag")[-2] %>% plot_grid()
Conclusions:
norin.lme2c %>% performance::check_model()
## Could not compute standard errors from random effects for diagnostic plot.
## Homogeneity of variance could not be computed. Cannot extract residual variance from objects of class 'lme'.
Unfortunately, this did not work…
## norin.resid = simulateResiduals(norin.lme2c, plot=TRUE)
plot_model(norin.lmer2c, type = "diag")[-2] %>% plot_grid()
Conclusions:
norin.lmer2c %>% performance::check_model()
Conclusions:
norin.resid <- norin.lmer2c %>% simulateResiduals(plot = TRUE)
Conclusions:
plot_model(norin.glmmTMB2c, type = "diag")[-2] %>% plot_grid()
Conclusions:
norin.glmmTMB2c %>% performance::check_model()
Conclusions:
norin.resid <- norin.glmmTMB2c %>% simulateResiduals(plot = TRUE)
Conclusions:
norin.lme2c %>% plot_model(type = "eff", terms = c("SMR_contr", "TRIAL"), show.data = TRUE)
We can also explore the fixed effects (as a caterpillar plot).
norin.lme2c %>% plot_model(type = "est")
And similarly, the individual random effects (intercepts of each fish).
# plot_model(norin.lme3a, type='re')
norin.lme2c %>%
allEffects() %>%
plot(multiline = TRUE, ci.style = "bands")
## Unfortunately, the use of scale() seems to prevent partial.residuals being calculated
## norin.lme2c %>% allEffects() %>% plot(multiline=TRUE, ci.style='bands', partial.residuals=TRUE)
norin.lme2c %>%
ggpredict(c("SMR_contr", "TRIAL")) %>%
plot(add.data = TRUE)
norin.lme2c %>%
ggemmeans(~ SMR_contr * TRIAL) %>%
plot(add.data = TRUE)
norin.lmer2c %>%
plot_model(type = "eff", show.data = TRUE) %>%
plot_grid()
norin.lmer2c %>% plot_model(type = "eff", terms = c("SMR_contr", "TRIAL"), show.data = TRUE)
We can also explore the fixed effects (as a caterpillar plot).
norin.lmer2c %>% plot_model(type = "est")
And similarly, the individual random effects (intercepts of each fish).
norin.lmer2c %>% plot_model(type = "re")
norin.lmer2c %>%
allEffects() %>%
plot(multiline = TRUE, ci.style = "bands")
norin.lmer2c %>%
ggpredict(c("SMR_contr", "TRIAL")) %>%
plot(add.data = TRUE)
norin.lmer2c %>%
ggemmeans(~ SMR_contr * TRIAL) %>%
plot(add.data = TRUE)
norin.glmmTMB2c %>%
plot_model(type = "eff", show.data = TRUE) %>%
plot_grid()
norin.glmmTMB2c %>%
plot_model(type = "eff", terms = c("SMR_contr", "TRIAL"), show.data = TRUE)
We can also explore the fixed effects (as a caterpillar plot).
norin.glmmTMB2c %>% plot_model(type = "est")
And similarly, the individual random effects (intercepts of each fish).
norin.glmmTMB2c %>% plot_model(type = "re")
norin.glmmTMB2c %>%
allEffects() %>%
plot(multiline = TRUE, ci.style = "bands")
norin.glmmTMB2c %>%
ggpredict(c("SMR_contr", "TRIAL")) %>%
plot(add.data = TRUE)
norin.glmmTMB2c %>%
ggemmeans(~ SMR_contr * TRIAL) %>%
plot(add.data = TRUE)
norin.lme2c %>% summary()
## Linear mixed-effects model fit by REML
## Data: norin
## AIC BIC logLik
## 1593.961 1635.028 -783.9803
##
## Random effects:
## Formula: ~TRIAL | FISHID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 16.323848 (Intr) TRIALH
## TRIALHypoxia 14.926864 -0.601
## TRIALLowSalinity 30.304028 -0.111 -0.215
## Residual 8.659992
##
## Fixed effects: CHANGE ~ TRIAL + scale(SMR_contr, scale = FALSE) + TRIAL:scale(SMR_contr, scale = FALSE)
## Value Std.Error DF
## (Intercept) 50.38244 2.385594 116
## TRIALHypoxia -50.22621 2.492663 116
## TRIALLowSalinity -42.22318 4.219647 116
## scale(SMR_contr, scale = FALSE) -21.55293 3.820913 58
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 13.51107 3.992402 116
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) 10.18925 6.758444 116
## t-value p-value
## (Intercept) 21.119453 0.0000
## TRIALHypoxia -20.149618 0.0000
## TRIALLowSalinity -10.006330 0.0000
## scale(SMR_contr, scale = FALSE) -5.640779 0.0000
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 3.384195 0.0010
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) 1.507633 0.1344
## Correlation:
## (Intr) TRIALH TRIALL s(Ss=F
## TRIALHypoxia -0.620
## TRIALLowSalinity -0.215 -0.035
## scale(SMR_contr, scale = FALSE) 0.000 0.000 0.000
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 0.000 0.000 0.000 -0.620
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) 0.000 0.000 0.000 -0.215
## TRIALHs=F
## TRIALHypoxia
## TRIALLowSalinity
## scale(SMR_contr, scale = FALSE)
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) -0.035
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.66585731 -0.30692682 0.02265975 0.29410064 1.72690964
##
## Number of Observations: 180
## Number of Groups: 60
Conclusions:
norin.lme2c %>% intervals()
## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) 45.657469 50.38244 55.10741
## TRIALHypoxia -55.163245 -50.22621 -45.28918
## TRIALLowSalinity -50.580718 -42.22318 -33.86563
## scale(SMR_contr, scale = FALSE) -29.201316 -21.55293 -13.90454
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 5.603613 13.51107 21.41852
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) -3.196697 10.18925 23.57520
##
## Random Effects:
## Level: FISHID
## lower est. upper
## sd((Intercept)) 5.0780215 16.3238483 52.4747730
## sd(TRIALHypoxia) 0.9252388 14.9268645 240.8149035
## sd(TRIALLowSalinity) 15.1035225 30.3040284 60.8026464
## cor((Intercept),TRIALHypoxia) -0.8169111 -0.6005325 -0.2359666
## cor((Intercept),TRIALLowSalinity) -0.8197871 -0.1109852 0.7321192
## cor(TRIALHypoxia,TRIALLowSalinity) -0.9832297 -0.2148081 0.9603211
##
## Within-group standard error:
## lower est. upper
## 0.1424309 8.6599921 526.5391714
norin.lme2c %>% tidy(conf.int = TRUE)
norin.lme2c %>%
tidy(conf.int = TRUE) %>%
kable()
| effect | group | term | estimate | std.error | df | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|---|
| fixed | fixed | (Intercept) | 50.3824384 | 2.385594 | 116 | 21.119453 | 0.0000000 | 45.6574691 | 55.1074076 |
| fixed | fixed | TRIALHypoxia | -50.2262117 | 2.492663 | 116 | -20.149618 | 0.0000000 | -55.1632454 | -45.2891780 |
| fixed | fixed | TRIALLowSalinity | -42.2231759 | 4.219647 | 116 | -10.006330 | 0.0000000 | -50.5807177 | -33.8656340 |
| fixed | fixed | scale(SMR_contr, scale = FALSE) | -21.5529273 | 3.820913 | 58 | -5.640779 | 0.0000005 | -29.2013157 | -13.9045389 |
| fixed | fixed | TRIALHypoxia:scale(SMR_contr, scale = FALSE) | 13.5110683 | 3.992402 | 116 | 3.384195 | 0.0009745 | 5.6036133 | 21.4185233 |
| fixed | fixed | TRIALLowSalinity:scale(SMR_contr, scale = FALSE) | 10.1892534 | 6.758444 | 116 | 1.507633 | 0.1343671 | -3.1966967 | 23.5752035 |
| ran_pars | FISHID | sd_(Intercept) | 16.3238483 | NA | NA | NA | NA | 5.0780215 | 52.4747730 |
| ran_pars | FISHID | cor_TRIALHypoxia.(Intercept) | -0.6005325 | NA | NA | NA | NA | -0.8169111 | -0.2359666 |
| ran_pars | FISHID | cor_TRIALLowSalinity.(Intercept) | -0.1109852 | NA | NA | NA | NA | -0.8197871 | 0.7321192 |
| ran_pars | FISHID | sd_TRIALHypoxia | 14.9268645 | NA | NA | NA | NA | 0.9252388 | 240.8149035 |
| ran_pars | FISHID | cor_TRIALLowSalinity.TRIALHypoxia | -0.2148081 | NA | NA | NA | NA | -0.9832297 | 0.9603211 |
| ran_pars | FISHID | sd_TRIALLowSalinity | 30.3040284 | NA | NA | NA | NA | 15.1035225 | 60.8026464 |
| ran_pars | Residual | sd_Observation | 8.6599921 | NA | NA | NA | NA | NA | NA |
Conclusions:
# warning this is only appropriate for html output
norin.lme2c %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)
| CHANGE | ||||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 50.38 | 2.39 | 45.66 – 55.11 | <0.001 |
| TRIAL [Hypoxia] | -50.23 | 2.49 | -55.16 – -45.29 | <0.001 |
| TRIAL [LowSalinity] | -42.22 | 4.22 | -50.58 – -33.87 | <0.001 |
| SMR contr | -21.55 | 3.82 | -29.20 – -13.90 | <0.001 |
|
TRIAL [Hypoxia] * SMR contr |
13.51 | 3.99 | 5.60 – 21.42 | 0.001 |
|
TRIAL [LowSalinity] * SMR contr |
10.19 | 6.76 | -3.20 – 23.58 | 0.134 |
| Random Effects | ||||
| σ2 | 75.00 | |||
| τ00 FISHID | 266.47 | |||
| τ11 FISHID.TRIALHypoxia | 222.81 | |||
| τ11 FISHID.TRIALLowSalinity | 918.33 | |||
| ρ01 | -0.60 | |||
| -0.11 | ||||
| ICC | 0.87 | |||
| N FISHID | 60 | |||
| Observations | 180 | |||
| Marginal R2 / Conditional R2 | 0.494 / 0.935 | |||
| AIC | 1593.961 | |||
norin.lmer2c %>% summary()
## Linear mixed model fit by REML ['lmerMod']
## Formula: CHANGE ~ TRIAL + scale(SMR_contr, scale = FALSE) + (TRIAL | FISHID) +
## TRIAL:scale(SMR_contr, scale = FALSE)
## Data: norin
## Control: lmerControl(check.nobs.vs.nRE = "ignore")
##
## REML criterion at convergence: 1568
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.30188 -0.42411 0.03131 0.40638 2.38624
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## FISHID (Intercept) 198.27 14.081
## TRIALHypoxia 86.41 9.296 -0.60
## TRIALLowSalinity 781.94 27.963 0.03 -0.64
## Residual 143.19 11.966
## Number of obs: 180, groups: FISHID, 60
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 50.382 2.386 21.119
## TRIALHypoxia -50.226 2.493 -20.150
## TRIALLowSalinity -42.223 4.220 -10.006
## scale(SMR_contr, scale = FALSE) -21.553 3.821 -5.641
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 13.511 3.992 3.384
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) 10.189 6.758 1.508
##
## Correlation of Fixed Effects:
## (Intr) TRIALH TRIALL s(Ss=F TRIALHs=F
## TRIALHypoxi -0.620
## TRIALLwSlnt -0.215 -0.035
## s(SMR_,s=FA 0.000 0.000 0.000
## TRIALH:(s=F 0.000 0.000 0.000 -0.620
## TRIALLS:s=F 0.000 0.000 0.000 -0.215 -0.035
Conclusions:
norin.lmer2c %>% confint()
## 2.5 % 97.5 %
## .sig01 0.000000 21.38850
## .sig02 NA NA
## .sig03 NA NA
## .sig04 NA NA
## .sig05 NA NA
## .sig06 NA NA
## .sigma NA NA
## (Intercept) 45.710775 55.05408
## TRIALHypoxia -55.107544 -45.34488
## TRIALLowSalinity -50.486428 -33.95992
## scale(SMR_contr, scale = FALSE) -29.035350 -14.07051
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 5.692828 21.32930
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) -3.045679 23.42418
norin.lmer2c %>% tidy(conf.int = TRUE)
norin.lmer2c %>%
tidy(conf.int = TRUE) %>%
kable()
| effect | group | term | estimate | std.error | statistic | conf.low | conf.high |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 50.3824384 | 2.385590 | 21.119486 | 45.706768 | 55.05811 |
| fixed | NA | TRIALHypoxia | -50.2262117 | 2.492656 | -20.149672 | -55.111729 | -45.34069 |
| fixed | NA | TRIALLowSalinity | -42.2231759 | 4.219658 | -10.006302 | -50.493554 | -33.95280 |
| fixed | NA | scale(SMR_contr, scale = FALSE) | -21.5529273 | 3.820907 | -5.640788 | -29.041768 | -14.06409 |
| fixed | NA | TRIALHypoxia:scale(SMR_contr, scale = FALSE) | 13.5110683 | 3.992391 | 3.384205 | 5.686126 | 21.33601 |
| fixed | NA | TRIALLowSalinity:scale(SMR_contr, scale = FALSE) | 10.1892534 | 6.758463 | 1.507629 | -3.057090 | 23.43560 |
| ran_pars | FISHID | sd__(Intercept) | 14.0807964 | NA | NA | NA | NA |
| ran_pars | FISHID | cor__(Intercept).TRIALHypoxia | -0.5968891 | NA | NA | NA | NA |
| ran_pars | FISHID | cor__(Intercept).TRIALLowSalinity | 0.0337670 | NA | NA | NA | NA |
| ran_pars | FISHID | sd__TRIALHypoxia | 9.2958585 | NA | NA | NA | NA |
| ran_pars | FISHID | cor__TRIALHypoxia.TRIALLowSalinity | -0.6361548 | NA | NA | NA | NA |
| ran_pars | FISHID | sd__TRIALLowSalinity | 27.9632565 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 11.9663529 | NA | NA | NA | NA |
Conclusions:
# warning this is only appropriate for html output
norin.lmer2c %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)
| CHANGE | ||||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 50.38 | 2.39 | 45.67 – 55.09 | <0.001 |
| TRIAL [Hypoxia] | -50.23 | 2.49 | -55.15 – -45.31 | <0.001 |
| TRIAL [LowSalinity] | -42.22 | 4.22 | -50.55 – -33.89 | <0.001 |
| SMR contr | -21.55 | 3.82 | -29.10 – -14.01 | <0.001 |
|
TRIAL [Hypoxia] * SMR contr |
13.51 | 3.99 | 5.63 – 21.39 | 0.001 |
|
TRIAL [LowSalinity] * SMR contr |
10.19 | 6.76 | -3.15 – 23.53 | 0.134 |
| Random Effects | ||||
| σ2 | 143.19 | |||
| τ00 FISHID | 198.27 | |||
| τ11 FISHID.TRIALHypoxia | 86.41 | |||
| τ11 FISHID.TRIALLowSalinity | 781.94 | |||
| ρ01 | -0.60 | |||
| 0.03 | ||||
| ICC | 0.76 | |||
| N FISHID | 60 | |||
| Observations | 180 | |||
| Marginal R2 / Conditional R2 | 0.494 / 0.877 | |||
| AIC | 1593.961 | |||
norin.glmmTMB2c %>% summary()
## Family: gaussian ( identity )
## Formula:
## CHANGE ~ TRIAL + scale(SMR_contr, scale = FALSE) + (TRIAL | FISHID) +
## TRIAL:scale(SMR_contr, scale = FALSE) + offset(scale(MASS,
## scale = FALSE))
## Data: norin
##
## AIC BIC logLik deviance df.resid
## 1606.3 1647.8 -790.1 1580.3 173
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## FISHID (Intercept) 369.934 19.234
## TRIALHypoxia 341.606 18.483 -0.54
## TRIALLowSalinity 1032.463 32.132 -0.17 -0.06
## Residual 8.349 2.889
## Number of obs: 180, groups: FISHID, 60
##
## Dispersion estimate for gaussian family (sigma^2): 8.35
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) 52.290 2.511 20.825
## TRIALHypoxia -56.118 2.444 -22.964
## TRIALLowSalinity -42.055 4.182 -10.057
## scale(SMR_contr, scale = FALSE) -19.840 4.022 -4.933
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 11.770 3.914 3.007
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) 10.884 6.698 1.625
## Pr(>|z|)
## (Intercept) < 2e-16 ***
## TRIALHypoxia < 2e-16 ***
## TRIALLowSalinity < 2e-16 ***
## scale(SMR_contr, scale = FALSE) 8.08e-07 ***
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 0.00264 **
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) 0.10414
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cov2cor(vcov(norin.glmmTMB2c)$cond)
## (Intercept) TRIALHypoxia
## (Intercept) 1.000000e+00 -5.472502e-01
## TRIALHypoxia -5.472502e-01 1.000000e+00
## TRIALLowSalinity -1.759538e-01 -3.991425e-02
## scale(SMR_contr, scale = FALSE) 3.764461e-15 -1.688551e-15
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) -9.655529e-16 -4.221495e-18
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) -8.500783e-16 6.952161e-16
## TRIALLowSalinity
## (Intercept) -1.759538e-01
## TRIALHypoxia -3.991425e-02
## TRIALLowSalinity 1.000000e+00
## scale(SMR_contr, scale = FALSE) -1.161186e-15
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 7.499722e-16
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) -2.198129e-15
## scale(SMR_contr, scale = FALSE)
## (Intercept) 2.642206e-15
## TRIALHypoxia -1.904170e-15
## TRIALLowSalinity -1.350437e-15
## scale(SMR_contr, scale = FALSE) 1.000000e+00
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) -5.472502e-01
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) -1.759538e-01
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)
## (Intercept) -1.892885e-15
## TRIALHypoxia 4.931878e-16
## TRIALLowSalinity 1.218256e-15
## scale(SMR_contr, scale = FALSE) -5.472502e-01
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 1.000000e+00
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) -3.991425e-02
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)
## (Intercept) -1.281649e-15
## TRIALHypoxia 9.201455e-16
## TRIALLowSalinity -3.060103e-15
## scale(SMR_contr, scale = FALSE) -1.759538e-01
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) -3.991425e-02
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) 1.000000e+00
norin.glmmTMB2c %>%
vcov() %>%
`[[`("cond") %>%
cov2cor() %>%
round(3)
## (Intercept) TRIALHypoxia
## (Intercept) 1.000 -0.547
## TRIALHypoxia -0.547 1.000
## TRIALLowSalinity -0.176 -0.040
## scale(SMR_contr, scale = FALSE) 0.000 0.000
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 0.000 0.000
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) 0.000 0.000
## TRIALLowSalinity
## (Intercept) -0.176
## TRIALHypoxia -0.040
## TRIALLowSalinity 1.000
## scale(SMR_contr, scale = FALSE) 0.000
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 0.000
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) 0.000
## scale(SMR_contr, scale = FALSE)
## (Intercept) 0.000
## TRIALHypoxia 0.000
## TRIALLowSalinity 0.000
## scale(SMR_contr, scale = FALSE) 1.000
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) -0.547
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) -0.176
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)
## (Intercept) 0.000
## TRIALHypoxia 0.000
## TRIALLowSalinity 0.000
## scale(SMR_contr, scale = FALSE) -0.547
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 1.000
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) -0.040
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)
## (Intercept) 0.000
## TRIALHypoxia 0.000
## TRIALLowSalinity 0.000
## scale(SMR_contr, scale = FALSE) -0.176
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) -0.040
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) 1.000
Conclusions:
norin.glmmTMB2c %>% confint()
## 2.5 % 97.5 %
## (Intercept) 4.736902e+01 5.721163e+01
## TRIALHypoxia -6.090730e+01 -5.132813e+01
## TRIALLowSalinity -5.025118e+01 -3.385951e+01
## scale(SMR_contr, scale = FALSE) -2.772227e+01 -1.195774e+01
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 4.099059e+00 1.944164e+01
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) -2.242694e+00 2.401121e+01
## Std.Dev.(Intercept) 1.783349e-05 2.074376e+07
## Std.Dev.TRIALHypoxia 1.589532e-12 2.149100e+14
## Std.Dev.TRIALLowSalinity 1.525316e-03 6.768844e+05
## Cor.TRIALHypoxia.(Intercept) -9.940613e-01 9.919646e-01
## Cor.TRIALLowSalinity.(Intercept) -6.845945e-01 6.843812e-01
## Cor.TRIALLowSalinity.TRIALHypoxia 6.014078e-01 7.708813e-01
## Estimate
## (Intercept) 52.2903273
## TRIALHypoxia -56.1177117
## TRIALLowSalinity -42.0553426
## scale(SMR_contr, scale = FALSE) -19.8400009
## TRIALHypoxia:scale(SMR_contr, scale = FALSE) 11.7703509
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) 10.8842586
## Std.Dev.(Intercept) 19.2336575
## Std.Dev.TRIALHypoxia 18.4825955
## Std.Dev.TRIALLowSalinity 32.1319577
## Cor.TRIALHypoxia.(Intercept) -0.5432691
## Cor.TRIALLowSalinity.(Intercept) -0.1658526
## Cor.TRIALLowSalinity.TRIALHypoxia -0.0552648
norin.glmmTMB2c %>% tidy()
## There seems to be a bug, doesn't work with random effects here
norin.glmmTMB2c %>% tidy(effects = "fixed", conf.int = TRUE)
norin.glmmTMB2c %>%
tidy(effects = "fixed", conf.int = TRUE) %>%
kable()
| effect | component | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|
| fixed | cond | (Intercept) | 52.29033 | 2.510916 | 20.825196 | 0.0000000 | 47.369022 | 57.21163 |
| fixed | cond | TRIALHypoxia | -56.11771 | 2.443711 | -22.964138 | 0.0000000 | -60.907297 | -51.32813 |
| fixed | cond | TRIALLowSalinity | -42.05534 | 4.181626 | -10.057175 | 0.0000000 | -50.251179 | -33.85951 |
| fixed | cond | scale(SMR_contr, scale = FALSE) | -19.84000 | 4.021637 | -4.933314 | 0.0000008 | -27.722265 | -11.95774 |
| fixed | cond | TRIALHypoxia:scale(SMR_contr, scale = FALSE) | 11.77035 | 3.913997 | 3.007246 | 0.0026363 | 4.099058 | 19.44164 |
| fixed | cond | TRIALLowSalinity:scale(SMR_contr, scale = FALSE) | 10.88426 | 6.697548 | 1.625111 | 0.1041389 | -2.242694 | 24.01121 |
Conclusions:
# warning this is only appropriate for html output
norin.glmmTMB2c %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)
| CHANGE | ||||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 52.29 | 2.51 | 47.37 – 57.21 | <0.001 |
| TRIAL [Hypoxia] | -56.12 | 2.44 | -60.91 – -51.33 | <0.001 |
| TRIAL [LowSalinity] | -42.06 | 4.18 | -50.25 – -33.86 | <0.001 |
| SMR contr | -19.84 | 4.02 | -27.72 – -11.96 | <0.001 |
|
TRIAL [Hypoxia] * SMR contr |
11.77 | 3.91 | 4.10 – 19.44 | 0.003 |
|
TRIAL [LowSalinity] * SMR contr |
10.88 | 6.70 | -2.24 – 24.01 | 0.104 |
| Random Effects | ||||
| σ2 | 8.35 | |||
| τ00 FISHID | 369.93 | |||
| τ11 FISHID.TRIALHypoxia | 341.61 | |||
| τ11 FISHID.TRIALLowSalinity | 1032.46 | |||
| ρ01 | -0.54 | |||
| -0.17 | ||||
| ICC | 0.99 | |||
| N FISHID | 60 | |||
| Observations | 180 | |||
| Marginal R2 / Conditional R2 | 0.501 / 0.993 | |||
| AIC | 1606.275 | |||
The presence of a significant interaction suggests that the nature of the relationship between SMR change and the control SMR differs between the three trials (High Temperature, Hypoxia and Low Salinity). Equivalently, the differences between the three trials is not constant over all levels of control SMR.
To further tease these interactions apart, we could either:
We will now explore each of the above, although typically we would only do one or the other (depending on which variable was of more interest).
norin.lme2c %>%
emtrends(~TRIAL, var = "SMR_contr") %>%
pairs()
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia -13.51 3.99 116 -3.384 0.0028
## HighTemperature - LowSalinity -10.19 6.76 116 -1.508 0.2911
## Hypoxia - LowSalinity 3.32 7.97 116 0.417 0.9088
##
## Degrees-of-freedom method: containment
## P value adjustment: tukey method for comparing a family of 3 estimates
Conclusions:
We will estimated the pairwise differences between each of the three trials at the minimum, mean and maximum control SMR.
norin.grid <- with(norin, list(SMR_contr = Hmisc::smean.sdl(SMR_contr)))
norin.grid
## $SMR_contr
## Mean Lower Upper
## 5.178857 3.926671 6.431044
norin.lme2c %>%
emmeans(~ TRIAL | SMR_contr, at = norin.grid) %>%
pairs()
## SMR_contr = 5.18:
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia 50.23 2.49 116 20.150 <.0001
## HighTemperature - LowSalinity 42.22 4.22 116 10.006 <.0001
## Hypoxia - LowSalinity -8.00 4.98 116 -1.608 0.2462
##
## SMR_contr = 3.93:
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia 67.14 5.59 116 12.020 <.0001
## HighTemperature - LowSalinity 54.98 9.46 116 5.814 <.0001
## Hypoxia - LowSalinity -12.16 11.15 116 -1.091 0.5217
##
## SMR_contr = 6.43:
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia 33.31 5.59 116 5.963 <.0001
## HighTemperature - LowSalinity 29.46 9.46 116 3.116 0.0065
## Hypoxia - LowSalinity -3.84 11.15 116 -0.345 0.9366
##
## Degrees-of-freedom method: containment
## P value adjustment: tukey method for comparing a family of 3 estimates
Conclusions:
norin.lme2c %>% r.squaredGLMM()
## R2m R2c
## [1,] 0.4942092 0.935456
## Nakagawa's R2
norin.lme2c %>% performance::r2_nakagawa()
## # R2 for Mixed Models
##
## Conditional R2: 0.935
## Marginal R2: 0.494
norin.lmer2c %>%
emtrends(~TRIAL, var = "SMR_contr") %>%
pairs()
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia -13.51 3.99 58 -3.384 0.0036
## HighTemperature - LowSalinity -10.19 6.76 58 -1.508 0.2949
## Hypoxia - LowSalinity 3.32 7.97 58 0.417 0.9088
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
norin.emt <- norin.lmer2c %>%
emtrends(~TRIAL, var = "SMR_contr") %>%
pairs() %>%
as.data.frame()
Conclusions:
norin.grid <- with(norin, list(SMR_contr = Hmisc::smean.sdl(SMR_contr)))
norin.grid
## $SMR_contr
## Mean Lower Upper
## 5.178857 3.926671 6.431044
norin.lmer2c %>%
emmeans(~ TRIAL | SMR_contr, at = norin.grid) %>%
pairs()
## SMR_contr = 5.18:
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia 50.23 2.49 58 20.150 <.0001
## HighTemperature - LowSalinity 42.22 4.22 58 10.006 <.0001
## Hypoxia - LowSalinity -8.00 4.98 58 -1.608 0.2503
##
## SMR_contr = 3.93:
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia 67.14 5.59 58 12.020 <.0001
## HighTemperature - LowSalinity 54.98 9.46 58 5.814 <.0001
## Hypoxia - LowSalinity -12.16 11.15 58 -1.091 0.5236
##
## SMR_contr = 6.43:
## contrast estimate SE df t.ratio p.value
## HighTemperature - Hypoxia 33.31 5.59 58 5.963 <.0001
## HighTemperature - LowSalinity 29.46 9.46 58 3.116 0.0079
## Hypoxia - LowSalinity -3.84 11.15 58 -0.345 0.9367
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
Conclusions:
norin.lmer2c %>% r.squaredGLMM()
## R2m R2c
## [1,] 0.4942089 0.876762
## Nakagawa's R2
norin.lmer2c %>% performance::r2_nakagawa()
## # R2 for Mixed Models
##
## Conditional R2: 0.877
## Marginal R2: 0.494
norin.glmmTMB2c %>%
emtrends(~TRIAL, var = "SMR_contr") %>%
pairs() %>%
summary(infer = TRUE)
norin.emt <- norin.glmmTMB2c %>%
emtrends(~TRIAL, var = "SMR_contr") %>%
pairs() %>%
as.data.frame()
Conclusions:
norin.grid <- with(norin, list(SMR_contr = Hmisc::smean.sdl(SMR_contr)))
norin.grid
## $SMR_contr
## Mean Lower Upper
## 5.178857 3.926671 6.431044
norin.glmmTMB2c %>%
emmeans(~ TRIAL | SMR_contr, at = norin.grid) %>%
pairs() %>%
summary(infer = TRUE)
Conclusions:
norin.glmmTMB2c %>% r.squaredGLMM()
## R2m R2c
## [1,] 0.501083 0.9934839
## Nakagawa's R2
norin.glmmTMB2c %>% performance::r2_nakagawa()
## # R2 for Mixed Models
##
## Conditional R2: 0.993
## Marginal R2: 0.501
norin.grid <- with(norin, list(SMR_contr = modelr::seq_range(SMR_contr, n = 100)))
newdata <- norin.lme2c %>%
emmeans(~ SMR_contr | TRIAL, at = norin.grid) %>%
as.data.frame()
head(newdata)
ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
geom_line(aes(, color = TRIAL)) +
theme_classic() +
theme(
legend.position = c(0.99, 0.99),
legend.justification = c(1, 1)
)
## The .fixed values are the predicted values without random effects
obs <- norin.lme2c %>%
augment() %>%
mutate(PartialObs = .fixed + .resid)
ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
geom_line(aes(, color = TRIAL)) +
geom_point(data = obs, aes(y = PartialObs, color = TRIAL)) +
geom_point(data = norin, aes(y = CHANGE), color = "gray") +
theme_classic()
g1 <- ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
geom_line(aes(, color = TRIAL)) +
geom_point(data = obs, aes(y = PartialObs, color = TRIAL)) +
theme_classic()
g2 <- ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
geom_line(aes(, color = TRIAL)) +
geom_point(data = norin, aes(y = CHANGE, color = TRIAL)) +
theme_classic()
g1 + g2
norin.grid <- with(norin, list(SMR_contr = modelr::seq_range(SMR_contr, n = 100)))
newdata <- norin.lmer2c %>%
emmeans(~ SMR_contr | TRIAL, at = norin.grid) %>%
as.data.frame()
head(newdata)
ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
geom_line(aes(, color = TRIAL)) +
theme_classic() +
theme(
legend.position = c(0.99, 0.99),
legend.justification = c(1, 1)
)
## The .fixed values are the predicted values without random effects
obs <- norin.lmer2c %>%
augment() %>% ## unfortunately, augment.lmer does not back-transform the scale()
dplyr::rename(SMR_contr = contains("SMR_contr")) %>%
mutate(
PartialObs = .fixed + .resid,
unscale(SMR_contr),
SMR_contr = V1
) # because ggfortify::unscale() returns a data.frame
ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
geom_line(aes(, color = TRIAL)) +
geom_point(data = obs, aes(y = PartialObs, color = TRIAL)) +
geom_point(data = norin, aes(y = CHANGE), color = "gray") +
theme_classic()
g1 <- ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
geom_line(aes(, color = TRIAL)) +
geom_point(data = obs, aes(y = PartialObs, color = TRIAL)) +
theme_classic()
g2 <- ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
geom_line(aes(, color = TRIAL)) +
geom_point(data = norin, aes(y = CHANGE, color = TRIAL)) +
theme_classic()
g1 + g2
norin.grid <- with(norin, list(SMR_contr = modelr::seq_range(SMR_contr, n = 100)))
newdata <- norin.glmmTMB2c %>%
emmeans(~ SMR_contr | TRIAL, at = norin.grid) %>%
as.data.frame()
head(newdata)
ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
geom_line(aes(, color = TRIAL)) +
theme_classic() +
theme(
legend.position = c(0.99, 0.99),
legend.justification = c(1, 1)
)
## The .fixed values are the predicted values without random effects
obs <- norin %>%
mutate(
.fixed = predict(norin.glmmTMB2c, re.form = NA),
.resid = residuals(norin.glmmTMB2c),
PartialObs = .fixed + .resid
)
ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
geom_line(aes(, color = TRIAL)) +
geom_point(data = obs, aes(y = PartialObs, color = TRIAL)) +
geom_point(data = norin, aes(y = CHANGE), color = "gray") +
theme_classic()
g1 <- ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
geom_line(aes(, color = TRIAL)) +
geom_point(data = obs, aes(y = PartialObs, color = TRIAL)) +
theme_classic()
g2 <- ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
geom_line(aes(, color = TRIAL)) +
geom_point(data = norin, aes(y = CHANGE, fill = TRIAL), color = "black", shape = 21) +
theme_classic()
g1 + g2